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Abstract 

The effect of artificial diffusion on discrete shock structures is examined for a family of schemes which includes 
scalar diffusion, convective upwind and split pressure (CUSP) schemes, and upwind schemes with characteristic 
splitting. The analysis leads to conditions on the diffusive flux such that stationary discrete shocks can contain 
a single interior point. The simplest formulation which meets these conditions is a CUSP scheme in which 
the coefficients of the pressure differences is fully determined by the coefficient of convective diffusion. It is 
also shown how both the characteristic and CUSP schemes can be modified to preserve constant stagnation 
enthalpy in steady flow, leading to four variants, the E and H-characteristic schemes, and the E and H-CUSP 
schemes. Numerical results are presented which confirm the properties of these schemes. 


1 Introduction 

The development of computational methods for the solution of gas dynamic equations has presented a contin- 
uing challenge. The goal of combining 

1. high accuracy 

2. high resolution of shock waves and contact discontinuities without oscillation 

3. minimum computational complexity 

in a single scheme has proved elusive. Two main issues in the design of non-oscillatory high resolution schemes 
were identified in the previous paper of this series [5]; first the design of scalar discrete schemes which guarantee 
the preservation of positivity and monotonicity in the solution, and second the construction of numerical fluxes 
for systems of equations to allow the proper resolution of complex wave interactions which may lead to the 
formation of both shock waves and contact discontinuities. The earlier paper develops systematic procedures 
for the design of scalar discretization schemes which satisfy positivity constraints. The present paper focuses 
on the design of numerical fluxes for the gas dynamic equations. 

Results presented in the previous paper confirm that stationary shocks can be resolved with a single interior 
point by combining either a symmetric limited positive (SLIP) scheme or an upstream limited positive (USLIP) 
scheme with a characteristic decomposition of the diffusive flux. The present paper presents an analysis of 
the conditions under which a discrete stationary shock can contain a single interior point. It emerges that 
a characteristic decomposition is not necessary to meet these conditions. Perfect single point discrete shocks 
completely free of oscillations can be produced by simpler flux splittings belonging to the class of convective 
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upwind and split pressure (CUSP) schemes, in which scalar diffusion is augmented by pressure differences. It 
is actually possible to obtain high resolution with almost no oscillation by introducing the right amount of of 
scalar diffusion, though this seems to result in a scheme which is less robust than the CUSP scheme. 

Section 2 reviews the shock jump conditions for one-dimensional flow, and their relationship to Roe’s lin- 
earization [9]. Section 3 reviews alternative splittings for a family of schemes. In all of them the diffusive flux 
is defined by a matrix which can be expressed as a polynomial function of the Jacobian matrix. Section 4 
examines semi-discrete schemes for the one-dimensional gas dynamic equations, and analyzes the conditions 
under which the numerical fluxes can be in perfect equilibrium when the discrete shock structure contains one 
interior point. These constraints can be satisfied by any numerical flux such that the equilibrium across the 
interface at the exit of the shock corresponds to the Hugoniot equation for a moving shock, while equilibrium 
across the interface at the entrance to the shock is maintained by full upwinding. The Roe linearization can 
be used to construct a variety of fluxes with these properties, with or without characteristic decomposition. In 
steady-state calculations the total enthalpy should be constant. Unfortunately numerical fluxes derived from 
the standard characteristic decomposition are not compatible with this property. Section 5 shows how the 
splittings can be modified so that this property is recovered while the discrete shock structure still has a single 
interior point. Section 6 discusses the implementation of limiters for these schemes. Numerical results which 
confirm the properties of the schemes are presented in section 7. 


2 Shock Jump Conditions and Roe Linearization 


The general one dimensional conservation law for a system of equations can be expressed as 

dw d 

Ht + = °' 

For the equations of gas dynamics the state and the flux vectors are 
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where p is the density, u is the velocity, E is the total energy, p is the pressure, and H is the stagnation 
enthalpy. If y is the ratio of specific heats and c is the speed of sound then 
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In a steady flow H is constant. This remains true for the discrete scheme only if the numerical diffusion is 
constructed so that it is compatible with this condition. 

When the flow is smooth it can be represented by the quasi-linear form 


dw dw 

~m +A ^ 


= 0 . 


where A(w) = , and the eigenvalues «, u+c and u — c of the Jacobian matrix A are the wave speeds for the 

three characteristics. In smooth flow the entropy S is constant along streamlines, and in isentropic flow the 
Riemann invariants JZ ± = u ± -~j are constant along the characteristics ^ = u ± c. These conditions may 
be expressed by the three equations, 
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Taking the dependent variables as 
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the equations can be expressed in the symmetric form 


dw 7 / ^ N dw 
— + yl(„)— = °, 


where 
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Depending on the initial data, there may not be a smooth solution of the conservation law (1). Nonlinear 
wave interactions along converging characteristics may lead to the formation and propagation of shock waves, 
while contact discontinuities may also appear. Denote the left and right states across a shock by subscripts L 
and R ) and let [/] and [w] be the jumps /r — fi and wr — w The shock jump condition is then 


[/] = s M , 


where S is the shock speed. 

In order to simplify the analysis of the equations when there are finite jumps in w and /, Roe introduced 
the linearization 

In ~ fi — w l )(w r — wi). 

where Arl{wr, wi) is a Jacobian matrix calculated from the left and right states in such a way that this 
relation is exact [9]. He showed that one way to do this is to introduce weighted averages 

u - + \fpL U L _ t/PrHr -h y/pZHL ^ 

\J~P~R d" \f~P~L ’ \/pR + y/pL 

into all the formulas in the standard expression for the Jacobian matrix A(w). In the case of a shock wave it 
now follows that 

Arl(wr — wi) — S(wr — wi). 

Thus the shock speed S is an eigenvalue of Arl , and the jump wr — wl is an eigenvector. In the case of a 
stationary shock 5 = 0. If we consider flow with u > 0 only the eigenvalue u — c can be zero. It follows that 
when u and c are calculated with Roe averages, u — c for a stationary shock. 


3 Alternative Splittings 


Suppose that the conservation law (1) is approximated over the interval (0, L) on a mesh with an interval Ax 
by the semi-discrete scheme 


a 

A ‘-dT + h i+i 


= 0 , 


( 3 ) 


where Wj denotes the value of the discrete solution in cell j , and hj + i is the numerical flux between cells j 
and j + 1. Let fj denote the flux vector f(wj ) evaluated for the state Wj . Suppose also that the numerical 
flux is 


h j+? - 2^ j+l + ^ ~ 


( 4 ) 


3 



where is a diffusive flux which is introduced to enable the scheme to resolve discontinuities without 

producing oscillations in the discrete solution. The diffusive flux is assumed to have the form 


d i+i = 2 a i+\ B i+\( w i+i- w i), 

where the matrix Sj+i determines the properties of the scheme, and the scaling factor is included for 

convenience. Introducing a Roe linearization, let + Wj) be an estimate of the Jacobian matrix 

with the property that 

A j+ ±(u>j+i - wj) = fj+i - fj. (5) 


+ ^ can be decomposed as 


A i+\ 


TAT- 1 , 


where the columns of T are the eigenvectors of Aj + ± , and A is a diagonal matrix containing its eigenvalues. 
Then the upwind scheme is produced by setting 


b j+\ = 


'i+i 


= T\\\T~ 


( 6 ) 


where the notation A^i 


is used to represent the absolute value of Aj+^ which is defined to be the matrix 
obtained by replacing the eigenvalues by their absolute values. Scalar diffusion is produced by setting 


B j+i = /• ( 7 ) 

An intermediate class of schemes can be formulated by defining the first order diffusive flux as a combination 
of differences of the state and flux vectors 


d j + i = a j+i (to i+ i - Wj) + 0 j+i (f j+1 - fj) . 

Schemes of this class are fully upwind in supersonic flow if one takes = 0 and ^ + i = sign (M) when the 
absolute value of the Mach number M exceeds 1. The flux vector / can be decomposed as 


/ = uw + fp , 


( 8 ) 


where 


Then 



( 9 ) 


fj + 1 “ fj = « K+l -Wj) + w (u j+ i - Uj) + f Pj+1 - f Pj , (10) 

where ti and w are the arithmetic averages 

- 1 / x 1 , 

u = 2 * Wi+1 + 2 + 

All these schemes can be obtained by representing + i as a polynomial in the matrix Aj + ± defined by 
equation (5). According to the Cayley-Hamilton theorem, a matrix satisfies its own characteristic equation. 
Therefore the third and higher powers of A can be eliminated, and there is no loss of generality in limiting 
Bj + 1 to a polynomial of degree 2, 

B j+i = a 0 I + a x A j+ i +a 2 ^| + i- (11) 

The characteristic upwind scheme for which + 1 = |a^ + i| is obtained by substituting Aj + 1 = TAT -1 , 
A 2 +1 = TA 2 T -1 . Then a 0 , Qfi, and a 2 are determined from the three equations 

« 0 -f 01 Ajt + = |Ajt | , k = 1,2,3. 


The same representation remains valid for three-dimensional flow because A ^+ 1 still has only three distinct 
eigenvalues u, u + c, u — c. 

Since Wj+\ — Wj approximates the diffusive flux introduces an error proportional to the mesh width, 

and both these schemes will be first order accurate unless compensating anti-diffusive terms are introduced. 
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Figure 1: Shock structure for single interior point. 


4 Conditions for a Stationary Shock 

The model of a discrete shock which will be examined is illustrated in figure (1). Suppose that wl and wr 
are left and right states which satisfy the jump conditions for a stationary shock, and that the corresponding 
fluxes are fi = /(^l) and /r = /(wr). Since the shock is stationary fi — Jr- The ideal discrete shock has 
constant states w L to the left and wr to the right, and a single point with an intermediate value w A . The 
intermediate value is needed to allow the discrete solution to correspond to a true solution in which the shock 
wave does not coincide with an interface between two mesh cells. According to equation (1) 

L rL rT 

w(T)dx = I w(0)dx+ / (/rb ~ fLB)dt, 

Jo Jo 

where Jib and /rr are the fluxes at the left and right boundaries. Assuming that the boundary conditions 
are compatible with a steady solution containing a stationary shock, the location x s of the shock is fixed by 
this equation, since 

/ w(T)dx ~ x 3 w l -F (L - x s )w R . 

Jo 

Similarly in the semi-discrete system 

Y^Wj(T) = ^2 w i(0)+ [ (fRB-fLB)dt. ( 12 ) 

. . Jo 

Thus w j (T) has a value which is determined by the initial and boundary conditions, and in general it is 
not possible for this value to be attained by a discrete solution without an intermediate point, because then 
the sum would be quantized, increasing by wr — wr whenever the shock location is shifted one cell to the 
right. 

Three diffusion models of varying complexity are examined in the following paragraphs to determine their 
ability to support the ideal shock structure containing a single interior point. These correspond to one, two 
or three terms in equation (11). 


4.1 Case 1 Scalar Diffusion 

The first model is simple scalar diffusion with 

d j+ i = ^a J+ £(u> j+1 -«>,). 

Consider the equilibrium in the cell immediately to the right of the shock. Using subscripts AR and RR to 
denote the values at the cell boundaries, the outgoing flux is 

fiRR - ^(/ft + /rt) - yORRiwR - Mr) -- /«. 
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while the incoming flux is 


>*AR = ^ (/« +/a)~ ~ w a)- 

For equilibrium these must be equal. It follows that 


Jr ~ Ja + <* A r(Wr - w A ) — 0. 

This is the Hugoniot condition for a shock moving to the left with a speed a A r. Introduce a Roe linearization 
with a mean Jacobian matrix A A r(w a , wr) such that 


Jr — Ja = A ar (w r - w A ). 


Then wr — w A is an eigenvector of A A r corresponding to the eigenvalue — a A R, The eigenvalues of A A r are 
ti, u + c and u — c. If we consider flow to the right with u > 0, and u < c, a solution with positive numerical 
diffusion is obtained by taking a AR = |u — c\. Then the intermediate value w A must lie on a Hugoniot curve 
defined by the right state wr. 

When the corresponding equilibrium is considered for a cell immediately to the left of a shock wave in a 
flow moving to the left, it is found that the diffusion coefficient should be \u + c|. Both cases can be satisfied 
by taking a = min(|u + c| , |u — c|). In the neighborhood of a stagnation point the accuracy can be improved 
by taking a proportional to u to prevent the numerical diffusion becoming undesirably large. This suggests 
the strategy of using a diffusion coefficient proportional to the smallest eigenvalue, or 


« j+ i = min | A t |, 

where A* are the eigenvalues u, u + c, and ti — c of A To prevent the scheme from admitting stationary 
expansion shocks which would violate the entropy condition, the diffusion coefficient may be redefined as 


a 


;+ \ “ 


(13) 


where 


Ajt = 


|A*| if|A fc |>c 

i (e + lAii!) if | At | < (, 


(14) 


and c is a positive threshold proportional to c. Recent work of Aiso [1] has established that in the scalar case 
this modification of the viscosity is sufficient to guarantee that the discrete solution will satisfy the entropy 
condition. The usual strategy in schemes using scalar diffusion has been to make the diffusion coefficient 
proportional to the maximum eigenvalue of the Jacobian matrix in order to make sure that the numerical 
viscosity for each characteristic variable is large enough to satisfy the positivity condition. Numerical tests with 
the alternative strategy of using the smallest eigenvalue confirm that very sharp discrete shocks are obtained, 
and that the scheme is robust with a viscosity threshold of the type defined by equation (14). 

To determine whether scalar diffusion can exactly support an ideal discrete shock it is also necessary to 
examine the equilibrium in the cell immediately before the shock. In this case the numerical fluxes are 


/ill = /l i 


and 

hLA = 7 j(fA + Jl ) - 2 a LA(WA ~ w l)- 

For equilibrium it is necessary that 


Ja~ Jl~ oc L a{w a - w L ) = 0, 

which is the Hugoniot condition for a shock moving to the right with a speed ai A . Introducing the Roe 
linearization, w A — wl must now be an eigenvector of Ala- The transition from L to A, however, is less than 
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the full jump for a stationary shock for which it is known that Roe averaging results in u = c. Thus it may 
be expected that u > c, and the choice a la — « — c ~ |u — c| could still allow the equilibrium condition to be 
satisfied. Then wa lies on a Hugoniot curve defined by the left state wl- 

The question now arises whether an intermediate state wa can be found that simultaneously lies on Hugoniot 
curves defined by the left and right states wl and wr, where these two states themselves satisfy the Hugoniot 
condition for a steady shock. It turns out that this is not possible. Let v = j be the specific volume. Then all 
possible shocks connecting wl and wr must satisfy the Hugoniot relation 


Prvr - plvl 



( PR + Pl)(vl - 


(15) 


This establishes a locus on ap-t diagram of a family of shocks as the shock speed is varied. The single point 
shock structure requires wa to lie on the Hugoniot curves defined by wl and wr. The curve defined from wl 
is 


Pa^a ~ Plvl — 



(Pa + Pl){vl -v A ) t 


(16) 


while the curve from wr is 

y ~ 1 

PrVr - Pava = 2 (PR + Pa)(v a - vr). (17) 

These intersect only when wa = vvr or wl . To prove this note that (15) can be written as 


Prvr — Plvl = &(prvl - Plvr ), (18) 

where a = Similarly (16) and (17) yield 

VL - CUV A Vr - OtV A 

Pa = PL = PR ■ 

V A - otv L V A - otVR 

Thus va satisfies a quadratic equation which may be written as 

(Prvr - Plvl)va - <*(pr - Pl){v 2 a + v l vr) «+■ oc 2 v a {prvl - Plvr) - 0 . 


Substituting from equation (18) 


va (prvl - Plvr) - {pr - Pl)(v 2 a + v L v R ) + v a (prvr - Plvl) = 0, 


or 

(PR ~ Pl){ v a - vr)(v a ~ v L ) = 0 . 

If p L ^ pr this has only the solutions va = vl or v A = vr. Therefore it is concluded that scalar diffusion 
cannot support a perfect discrete shock with a single interior point. Calculations of one-dimensional flows 
reveal an oscillation of very small amplitude upstream of the shock. In multidimensional flows, however, these 
oscillations are essentially imperceptible. 


4.2 Case 2 Characteristic Upwind Scheme 


The second case to be examined is the upwind scheme which results from characteristic decomposition, with 


B i +* = 




. This case has been studied by Roe [10], and it is known that the upwind scheme admits ideal 


shocks. Assuming flow to the right with u > 0, the fluxes in the cell to the right of the shock are now 

hRR - Jr, 


and 

h A R ~ 2 (/* + / 'a ) — 2 \Aar\ ( w r - u;^), 

yielding equilibrium if 


{A ar - - wa) = T { A - |A|)T \w R - w A ) - 0. 
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With u < c this is satisfied by the negative eigenvalue u — c, and since wr — wa is the corresponding eigenvector, 
the Hugoniot equation 

A - Sa = S(w H - w A ) 

is satisfied for the shock speed S — u — c. Thus wa again lies on a Hugoniot curve. At the entrance to the 
shock the transition from wi to wa is less than the full transition from wl to wr for which u = c. Thus 
a structure is admitted in which u > c in the transition from L to A, with the consequence that the flux is 
calculated from the upwind state 

/ LA - 2 (/>» + A) - 2 Ala{u>l - tv a) = A 

and equilibrium is maintained. 

4.3 Case 3 Convective Upwind and Split Pressure (CUSP) Scheme 

Characteristic decomposition allows equilibrium to be established through full upwinding of the flux entering 
the transition layer, while the flux leaving the transition layer satisfies a Hugoniot equation. This can also be 
accomplished by a less complex scheme. Suppose that the diffusive flux is defined as 

d j+ i = ^a‘c(u ) j+1 - Wj) + /?(/> + ! - fj), 

where the factor c is included so that a* is dimensionless. Let M be the Mach number If the flow is 

c 

supersonic an upwind scheme is obtained by setting 

a* = 0, ft = sign(M). 

Introducing the Roe linearization, the Mach number is calculated from tz and c, and at the entrance to the 
shock a transition to an intermediate value wa is admitted with u > c and 

Sla = f*. + h) - 2 (/ a ~ A) = A- 

The fluxes leaving and entering the cell immediately to the right of the shock are now 

Irr — /* > 

and 

/ar = 2 (A + A) - -a*c(w R - w A ) - -/?(/« - f A ). 

These are in equilibrium if 

cx* c 

A - A + Y~ p WR ~ Wa ^ = 

This is the Hugoniot equation for a shock moving to the left with a speed Also, introducing the Roe 

linearization, 

Q* C 

(Ar A + ~ w a) = 0 . 

Thus wr — wa is an eigenvector of Ara and is the corresponding eigenvalue. Since the eigenvalues are 

tt, tt + c and tz — c, the only choice which leads to positive diffusion when u > 0 is tz — c, yielding the relationship 

Gf*c = (1 + fi){c - u), 0 < u < c. 

Thus ft is uniquely determined once a* is chosen, leading to a one-parameter family of schemes. The choice 
ft = M corresponds to the Harten-Lax-Van Leer (HLL) scheme [4, 3], which is extremely diffusive. 

The term ft(fn — /a) contributes to the diffusion of the convective terms. Suppose that the convective terms 
are separated by splitting the flux according to equations (8), (9) and (10). Then the total effective coefficient 
of convective diffusion is 

otc — a*c + ftu. 
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The choice ac = u leads to low diffusion near a stagnation point, and also leads to a smooth continuation 
of convective diffusion across the sonic line since o* = 0 and j3 — 1 when \M\ > 1. The scheme must also 
be formulated so that the cases of u > 0 and u < 0 are treated symmetrically. Using the notation M — 

X ± = u zb c, this leads to the diffusion coefficients 


a = \M\ 


(19) 

+ max(0,^) 

if 0 < M < 1 


-max(0.a±g)) 

if - 1 < M < 0 

(20) 

sign(M) 

if \M\ > 1. 



Near a stagnation point a may be modified to a = ^ if |M| is smaller than a threshold e. 

4.3.1 Criteria for a single point shock 

The analysis of these three cases shows that a discrete shock structure with a single interior point is supported 
by artificial diffusion that satisfies the two conditions that 

1. it produces an upwind flux if the flow is determined to be supersonic through the interface 

2. it satisfies a generalized eigenvalue problem for the exit from the shock of the form 

{Aar — olarBar) (wr — ™a) = 0 , 

where Aar is the linearized Jacobian matrix and Bar is the matrix defining the diffusion for the interface 
AR. These two conditions are satisfied by both the characteristic and CUSP schemes. Scalar diffusion does 
not satisfy the first condition. 


5 Schemes Admitting Constant Total Enthalpy in Steady Flow 


In steady flow the stagnation enthalpy H is constant, corresponding to the fact that the energy and mass 
equations are consistent when the constant factor H is removed from the energy equation. Discrete and semi- 
discrete schemes do not necessarily satisfy this property. In the case of a semi-discrete scheme expressed in 
viscosity form, equations (3) and (4), a solution with constant H is admitted if the viscosity for the energy 
equation reduces to the viscosity for the continuity equation with p replaced by pH . When the standard 
characteristic decomposition (6) is used, the viscous fluxes for p and pH which result from composition of the 
fluxes for the characteristic variables do not have this property, and H is not constant in the discrete solution. 
In practice there is an excursion of H in the discrete shock structure which represents a local heat source. In 
very high speed flows the corresponding error in the temperature may lead to a wrong prediction of associated 
effects such as chemical reactions. 

The source of the error in the stagnation enthalpy is the discrepancy between the convective terms 


p 

pu 

pH 


in the flux vector, which contain pH , and the state vector which contains pE . 
introducing a modified state vector 


w h ~ 


p 

pu 

pH 


This may be remedied by 
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Then one introduces the linearization 


/r - ?L = A h (w hR - w hL ). 

Here Ah may be calculated in the same way as the standard Roe linearization. On introducing the vector 

( yfP > 

V = Jpu 

VP” 

all quantities in both / and w k are products of the form VjV k which have the property that a finite difference 
A (vjVk) between left and right states can be expressed as 

A (vjv k ) = VjAv k + v k Avj 

where Vj is the arithmetic mean \{yjn + vji). Therefore, 

Aw = BAv, Af — CAv = CB~ l Aw, 

where B and C can be expressed in terms of appropriate mean values of the quantities Vj . 

Define 

_ yfpRUR + yfPL^L ^ 

yJ~PR + \J~PL ' y/PR + y/PL 

and 


c=;/( 7 -!)(//- y)- 


Then 


( 0 


Ah = 


V -uH 

The eigenvalues of Ah are u, A + and A - where 



A* = 


7+1 

2 7 


7+1 


u± v (_ 2 r u)2+ 


( 21 ) 


Note that A + and A have the same sign as «+c and u — c, and change sign at the sonic line u = ±c. The 
corresponding eigenvectors of Ah are the columns of 


T — 


1 1 1 
u A + A" 
4 H H 


Also the left eigenvectors of Ah are the rows of 


( (A + -A ~)H 


t -' = d 


0 


(A+-A-) 


-(xtH- \~\) H-\ -(A- -«) | , 

V -(*-£) A+-U 


where 


Then 


D = ( A+ - A ~){H - —) = (A+ - A-)—. 


Ah = TAT -1 , 


10 



where 


The same development can be carried out for multidimensional flows. For convenience the formulas for the 
general multi-dimensional case are presented in the Appendix. 

Using the modified linearization both the characteristic upwind scheme and the CUSP scheme can be 
reformulated as follows to admit steady solutions with constant H . 

5.1 Case 1 Characteristic Upwind Scheme 

The diffusion for the characteristic upwind scheme is now defined to be 

dj + 1 = \ A hj+ ^ (w j+i -wj), 


in which \Ah I is defined to be 


where T is the eigenvector matrix of Ah> and 


\A h \ = TAT~\ 


In order to show that the scheme admits a solution with constant H y split the diffusion into two parts 

a _ jO) ,j(2) 

where is the contribution from |u|, and where is the contribution from |A+| and |A“| . Then 

2 J *T 2 


(t?2i = \ U X+ X- 

H H 


1 1 1 


T-'Awh 


0 |A + | |A-| \ 

0 |A+|A+ |A“ | A- 

0 |A+| H \X~\H ) 


and the third element of equals the first element multiplied by H . Also 


% - 


(7 - 1) l«l 


1 1 1 \ / HAp-A(pH) 

u X+ A- 0 

4 H H ) \ 0 


(y-l)p\u\AH 


and this is zero if H is constant. Thus both contributions are consistent with a steady solution in which H 
is constant. The two variations of the characteristic splitting can conveniently be distinguished as the E and 
H-characteristic schemes. 
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5.2 Case 2 CUSP Scheme 


The diffusive flux is now expressed as 

d j+t = \<**cAw h + 

where A denotes the difference from j -f 1 to j. Again equilibrium at the entrance is established by upwinding, 
while equilibrium at the exit requires 

. . o*c . . » a*c . 

A / + TjTp Awh = ( Ah + T^Tp^ Awh = °- 

Therefore, —f+j; must be an eigenvalue of Ah, and in the case u > 0, positive diffusion is obtained by taking 

a*c - -(1 + P)\~ . 


Now the split is redefined as 


/ - uw h + fp, 


where 


and the diffusive flux can be expressed as 


fp — 


( ° 

V 

0 


d ; +i = -arc A Wh + 0WhAu + f3 A f p . 


Then a and (3 are defined as before by equations (19) and (20), using the modified eigenvalues X ± defined 
equation (21). This splitting corresponds to the Liou-Steffen splitting [7, 13]. The splitting in which the 
convective terms contain pE corresponds to the wave particle splitting [8, 2]. As in the case of characteristic 
splitting, the two variations can conveniently be distinguished as the E-CUSP and H-CUSP schemes. 


6 Implementation of limiters for the CUSP scheme 


In the case of a scalar conservation law, high resolution schemes which guarantee the preservation of the 
positivity or monotonicity of the solution can be constructed by limiting the action of higher order or anti- 
diffusive terms, which might otherwise cause extrema to grow. Typically, these schemes, such as both the 
symmetric and upstream limited positive (SLIP and USLIP) schemes discussed in the previous paper in this 
series [5], compare the slope of the solution at nearby mesh intervals. The characteristic upwind scheme 
essentially applies the same construction to the characteristic variables, so that the solution is subject to 
controls on the formulation or growth of extrema of these variables. The fluxes appearing in the CUSP scheme 
have different slopes approaching from either side of the sonic line, and use of limiters which depends on 
comparisons of the slopes of these fluxes can lead to a loss of smoothness in the solution at the entrance to 
supersonic zones in the flow. 

An alternative formulation which avoids this difficulty, and may be used with either the characteristic 
upwind or the CUSP scheme, is to form the diffusive flux from left and right states at the cell interface. These 
are interpolated or extrapolated from nearby data, subject to limiters to preserve monotonicity, in a similar 
manner to the reconstruction of the solution in Van Leer’s MUSCL scheme [6]. Let 


R(u, u) = 1 


u — v 


M + HI 

where q is a positive power. Then R(u,v) = 0 when it and v have opposite' sign. Also define 

+ v). 


( 22 ) 


(23) 
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Let denote the kth element of the state vector w . Now define left and right states for each dependent 
variable separately as 


w 


(*) _ 


= w 


(k) + -HAw^,Aw^) 


„(*) _ 


= w® - \ L ( Aw f+§ > Au, !-V’ 


where 


A w j+ i = tfj + i 


Then 


W-W- 


= Auf ’, -L(A^,A W ^_ ) 1 ) 

J+2 J ^ 2 J 2 

and in the case of a scalar equation the SLIP scheme [5] is recovered by making the diffusive flux proportional 
to this difference. To implement the CUSP scheme the pressures p L and p R for the left and right states are 
determined from wl and w R . Then the diffusive flux is calculated by substituting wl for Wj and wr for + i 
to give 

1 1 

d j+ x = ~ot*c(w R - w L ) + -/?(/( w R ) - f(w L )). 


Similarly the characteristic upwind scheme is implemented by calculating + ^ from wr and wi. An alter- 
native reconstruction is to set 


w 

w 

It has been found that essentially 
the two interpolation formulas. 


(*) 

L 

w 

R 


= + R(Aw[ k +s , Au/*^)Ati;^\ 

= w {k) - R(Aw {k \ , Auf k \ )Aw (k ^ . 

J JT 2 J 2 J ' 2 


similar results are obtained in numerical calculations of steady flows using 


7 Numerical Results 

Extensive numerical tests have been performed with the E and H-characteristic and the E and H-CUSP schemes 
to verify their properties. Results for one, two and three-dimensional flows are presented in sections 7.1, 7.2 
and 7.3. 

7.1 One dimensional shock 

In order to verify the discrete structure of stationary shocks with the various schemes, calculations were 
performed for a one dimensional problem with initial data containing left and right states compatible with the 
Rankine-Hugoniot conditions. An intermediate state consisting of the arithmetic average of the left and right 
states was introduced at a single cell in the center of the domain. With this intermediate state the system 
is not in equilibrium, and the time dependent equations were solved to find an equilibrium solution with a 
stationary shock wave separating the left and right states. Tables 1 through 4 shows the results for a shock 
wave at Mach 20 for the E-characteristic, H-characteristic, E-CUSP and H-CUSP schemes. In all cases the 
SLIP construction was used with the limiter defined by equations (22) and (23), and q — 3. The tables show 
the values of p, tt, //, p, M and the entropy S = log ^ - log j . All four schemes display a perfect one 

point shock structure. The entropy is zero to 4 decimal places upstream of the shock, and is constant to 4 
decimal places downstream of the shock. There is a slight excursion of the entropy at the interior point in the 
results for the H-characteristic and H-CUSP schemes. Correspondingly there is an excursion in the stagnation 
enthalpy at the interior point in the results for the E-characteristic and E-CUSP schemes. It may be noted 
that the mass, momentum and energy of the initial data are not compatible with the final equilibrium state. 
According to equation (12) the total mass, momentum and energy must remain constant if the outflow flux f R 
remains equal to the inflow flux f L . Therefore f R must be allowed to vary according to an appropriate outflow 
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boundary condition to allow the total mass, momentum and energy to be adjusted to values compatible with 
equilibrium. 

7.2 Airfoil calculations 

The results of transonic flow calculations using the H-characteristic and H-CUSP schemes are compared in 
figures (2-8). The E-characteristic and E-CUSP schemes produce results which are very similar to the results of 
the H-characteristic and H-CUSP schemes, with small deviations in stagnation enthalpy. These are eliminated 
by the H-characteristic and H-CUSP schemes. The limiter defined by equations (22) and (23) was again used 
with q = 3 in both schemes to define left and right states in the manner described in section 6. The H-CUSP 
scheme was simplified by replacing the Roe averages (2) by arithmetic averages, and using A* = u ± c in the 
formula (20) for /?. It was also found that the term tends to reduce the rate of convergence to a steady 

state. Therefore it was attenuated by the factor Jp*\ _ ps | ^ where ps is the pressure at sonic speed, and 
Pl and pr are the pressures to the left and right. When the flow crosses the sonic line ps lies between pl 
and pRj and this factor becomes unity. Thus the full scheme is restored at a shock wave. All the calculations 
were performed with the five stage modified Runge-Kutta time stepping scheme described in reference [5]. 
Convergence to a steady state was accelerated by the multigrid also described in reference [5], using W-cycles 
in which a single time step is performed on each grid level during the descent towards coarser grids. The total 
amount of work in each W-cycle is about the same as two time steps on the fine grid. 

Calculations are presented for two well known airfoils, the RAE 2822 and the NACA 0012. The equations 
were discretized on meshes with O-topology extending out to a radius of about 100 chords. In each case the 
calculations were performed on a sequence of successively finer meshes from 40x8 to 320x64 cells, while the 
multigrid cycles on each of these meshes descended to a coarsest mesh of 10x2 cells. Figure 2 shows the inner 
parts of the 160x32 meshes for the two airfoils. Figures 3-8 show the final results for each scheme on 320x64 
meshes for the RAE 2822 airfoil at Mach .75 and 3° angle of attack, and for the NACA 0012 airfoil at Mach .8 
and 1.25° angle of attack, and also at Mach .85 and 1° angle of attack. In each case the convergence history is 
shown for 100 or 200 cycles, while the pressure distribution is displayed after a sufficient number of cycles for 
its convergence. The pressure distribution of the RAE 2822 airfoil converged in only 25 cycles. Convergence 
was slower for the NACA 0012 airfoil. In the case of flow at Mach .8 and 1.25° angle of attack, additional 
cycles were needed to damp out a wave downstream of the weak shock wave on the lower surface. 

As a further check on accuracy the drag coefficient should be zero in subsonic flow, or in shock free transonic 
flow. Tables 5 and 6 show the computed drag coefficient with the H-characteristic and H-CUSP schemes on a 
sequence of three meshes for three examples. The first two are subsonic flows over the RAE 2822 and NACA 
0012 airfoils at Mach .5 and 3° angle of attack. The third is the flow over the shock free Korn airfoil at its 
design point of Mach .75 and 0° angle of attack. The computed drag coefficients are slightly lower with the 
H-CUSP scheme: in all three cases the drag coefficient is calculated to be zero to four digits on a 160x32 mesh. 


7.3 Three-dimensional calculations for a swept wing 

As a further test of the performance of the H-CUSP scheme, the flow past the ONERA M6 wing was calculated 
on a mesh with C-H topology and 192x32x48 = 294912 cells. Figure 9 shows the result at Mach .84 and 3.06° 
angle of attack. This again verifies the non-oscillatory character of the solution, and the sharp resolution of 
shock waves. In this case 50 cycles were sufficient for convergence of the pressure distributions. 


8 Conclusion 

It was shown in the first paper in this series [5] that the concept of local extremum diminishing (LED) 
schemes provides a convenient framework for the formulation of non-oscillatory shock capturing schemes for 
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I 

P 

u 

H 

P 


S 

1 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

2 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

3 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

4 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

5 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

6 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

7 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

8 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

9 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

10 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

11 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

12 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

13 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

14 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

15 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

16 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

17 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

18 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

19 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

20 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

21 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

22 

4.4232 

6.9992 

268.9099 

308.8882 

0.7079 

37.5269 

23 

5.9259 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

24 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

25 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

26 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

27 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

28 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

29 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

30 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

31 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

32 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

33 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

34 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

35 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

36 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

37 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

38 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

39 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

40 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

41 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 

42 

5.9260 

3.9930 

283.5064 

466.5228 

0.3803 

37.6369 


Table 1: Shock Wave at Mach 20: ^characteristic scheme 





I 

P 

u 

H 

P 

M 

S 

1 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 


1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 


1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 


1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 


1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

6 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

7 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

8 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

9 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

10 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

11 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

12 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

13 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

14 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

15 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

16 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

17 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

18 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

19 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

m \ 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

H 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

22 

4.2476 

7.1962 

283.5073 

312.6405 

0.7089 

40.2710 

23 

5.9259 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

24 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

25 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

26 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

27 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

28 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

29 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

30 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

31 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

32 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

33 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

34 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

35 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

36 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

37 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

38 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

39 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

40 

5.9260 

3.9930 

283.5073 

466.5208 

0.3803 

37.6371 

41 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 

42 

5.9260 

3.9930 

283.5074 

466.5208 

0.3803 

37.6371 


Table 2: Shock Wave at Mach 20: H-characteristic scheme 


16 








I 

P 

u 

H 

P 

M 

S 

1 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

2 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

3 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

4 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

5 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

6 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

7 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

8 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

9 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

10 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

11 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

12 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

13 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

14 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

15 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

16 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

17 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

18 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

19 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

20 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

21 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

22 

4.4011 

7.0503 

268.7344 

306.6689 

0.7138 

37.5200 

23 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

24 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

25 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

26 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

27 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

28 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

29 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

30 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

31 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

32 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

33 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

34 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

35 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

36 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

37 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

38 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

39 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

40 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

41 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 

42 

5.9259 

3.9935 

283.4970 

466.4894 

0.3804 

37.6357 


Table 3 : Shock Wave at Mach 20 : E-CUSP scheme 
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I 

P 

u 

H 

P 

M 

s 

1 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

2 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

3 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

4 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

5 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

6 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

7 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

8 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

9 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

10 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

11 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

12 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

13 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

14 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

15 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

16 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

17 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

18 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

19 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 


1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

H 

1.0000 

23.6643 

283.5000 

1.0000 

20.0000 

0.0000 

22 

4.1924 

7.3248 

283.4960 

307.4467 

0.7229 

40.3353 

23 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

24 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

25 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

26 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

27 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

28 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

29 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

30 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

31 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

32 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

33 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

34 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

35 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

36 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

37 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

38 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

39 

5.9259 

3.9935 

283.4960 

466.4889 

0.3804 

37.6355 

40 

5.9259 

3.9935 

283.4961 

466.4889 

0.3804 

37.6355 

41 

5.9259 

3.9935 

283.4961 

466.4889 

0.3804 

37.6355 

42 

5.9259 

3.9935 

283.4961 

466.4889 

0.3804 

37.6355 


Table 4: Shock Wave at Mach 20: H-CUSP scheme 


Mesh 

RAE 2822 
Mach .50 a 3° 

NACA 0012 
Mach .50 a 3° 

Korn Airfoil 
Mach .75 a 0° 


■ 


.0126 

80x16 


.0017 

.0026 

160x32 

.0002 

.0002 

.0001 


Table 5: Drag Coefficient on a sequence of meshes: H-characteristic scheme 
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Mesh 

RAE 2822 
Mach .50 a 3° 

NACA 0012 
Mach .50 a 3° 

Korn Airfoil 
Mach .75 a 0° 









160x32 





Table 6: Drag Coefficient on a sequence of meshes: H-CUSP scheme 


compressible flow calculations. In the case of scalar conservation laws the LED property can be secured by 
corresponding symmetric and upstream limited positive (SLIP and USLIP) schemes. 

The different scalar constructions can be combined with alternate numerical fluxes to provide a matrix of 
schemes for the gas dynamic equations. The property of supporting stationary discrete shocks with a single 
interior point is shared by the characteristic and CUSP schemes. Each of these schemes can be modified to 
preserve constant stagnation enthalpy in steady flows, giving four variants, the E and H-characteristic schemes, 
and the E and H-CUSP schemes. The CUSP schemes are inexpensive. They introduce a minimum amount of 
numerical diffusion as the Mach number approaches zero. They are therefore also appropriate for viscous flow 
calculations in which it is important not to contaminate the boundary layer. 

The theoretical properties of these schemes are verified by numerical calculations of one-dimensional, two- 
dimensional and three-dimensional flows. References [11] and [12] evaluate the accuracy and efficiency of some 
of these schemes for the calculation of viscous flows. 
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Appendix: Eigenvalues and Eigenvectors for Gas Dynamic Equa- 
tions 

The Euler equations which describe the three-dimensional flow of an inviscid gas can be written as 

4 / WdV+ f (F x dS x + F y dS y + F z dS z ) = 0 
dt Jet JdCl 


(24) 


where W is the state vector, F x , F y and F z are the flux vectors, and dS X} dS y and dS z are the projections of 
the surface element in the x, y and z coordinate directions. Let u, v and w be the velocity components and p, 
p, E and H the density, pressure, total energy and total enthalpy. Then 


W = < 


> 

p 


r y 

pu 


' 1 

pv 


s y 

pw 

pu 


puu -1- p 


pvu 


pwu 

i pv 

► , F x = < 

puv 

> , Fy = i 

pvv -h p 


pwv 

pw 


puw 


pvw 


pww + p 

. pE . 


„ puH j 


. P VH > 


. p wH , 


Also, 






where q is the flow speed, and c is the speed of sound, 


q* = U 2 + V 2 + w\ C 2 =^ 

P 

When flow is smooth it can be represented by the quasi-linear differential equation 

dW 


dW 

dt 


dW 

dx 


+ A x + A z -^— - 0 


dW 

dz 


where A x , A y and A z are the Jacobian matrices 


A -M* a a 

r ~ dW' y dW' dW 

Under a change of variables to a new state vector W , equation (28) is transformed to 

dW - dW - dW ~ dW n 

~dr + A *~fc +Ay -di + A *-d 7 = 0 


(25) 


(26) 


(27) 


(28) 
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where 


A x = MA X M 


-l 


Ay = MAyM~ 


A 2 = M A Z M~ 


and 


M = 


dW 

dw 


The finite volume discretization requires the evaluation of the flux through a face with vector area S, 

F = F X S X 4- FySy 4- F Z S Z . 

The corresponding Jacobian matrix is 


A = — S X A X 4- SyAy H~ 


Let S be the magnitude of the face area 


s 2 = sl 4- s; 4- si 


and n x , n y and n z be the components of the unit normal 

5, 


S X Uy 

n x — “jT ^ n y = g ) ^"2 — g 


Also let Q be the flux rate 


Q = uS x + 4- 1 oS z 

Then the Jacobian matrix can be decomposed as 

A = MAM" 1 

where 


M = 



1 

0 

0 

0 

0 

\ 


( 1 

0 

0 

0 


u 

1 

0 

0 

0 



—u 

1 

0 

0 


V 

0 

1 

0 

0 


, M~ l = 

— V 

0 

1 

0 


w 

0 

0 

1 

0 



—w 

0 

0 

1 


al 

2 

u 

V 

w 

i 

<Y-1 

J 


(7-1)4 

-(7 - 1)« 

1 

1 

-(7 - 


and 


A = 


/ <3 
o 
o 
0 


5, 

Q 

o 

0 


0 

Q 

o 


s* 

0 

0 

Q 


Also the further transformation 
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produces the symmetric form 


S x c 

A= | SyC 
S z c 

\ ° 

Then A, A, and A can be decomposed as 

A = RAR~\ 
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where the diagonal matix A contains the eigenvalues 
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Also the right eigenvectors of A, A, and A are the columns of 

R=DN , R = P~ 1 N, R = MDN 
and the left eigenvectors of A, A, and A are the rows of 

R~ l = AT'zr 1 , R~ l = N-'P, R~ l = N-'D-'M- 1 

where 
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The decomposition to A corresponds to the introduction of primitive variables, scaled by a diagonal matrix, 


dW = DdW, W = 


The decomposition to ,4 corresponds to the introduction of symmetrizing variables dW, defined in differential 
form, scaled by another diagonal matrix 

\ 



P 

\ 



l 

0 

0 

0 

0 

\ 


u 




0 

p 

0 

0 

0 



V 


b> 

ii 


0 

0 

p 

0 

0 



w 




0 

0 

0 

p 

0 



p 

) 



0 

0 

0 

0 

p 

) 


/ 


dW = DdW, W = 


d£ 

pc 

du 


dv 
dw 

^ dp — c 2 dp J 


D = 


i \ 
0 
0 
0 

\ 0 


0 

£ 

c 

0 

0 

0 


0 

0 

£ 

c 

0 

0 


0 

0 

0 

£ 

c 

0 


0 \ 
0 
0 
0 

-7* / 


Multiplying out MDN , the right eigenvectors of A can be expressed as follows. Set 
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Then the eigenvectors corresponding to the eigenvalue Q are 


r i = n x v o + ui, r 2 = n y v 0 + t> 2 , r 3 = n z v 0 + u 3 , 

Here vo represents an entropy wave and represent vorticity waves. Also let q n denote the normal velocity 
Q/S , and set 


1 ) 


0 1 
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\ In J 


Then the last two eigenvectors, corresponding to pressure waves, are 




V 4 -f cu 5 , r 5 = V 4 - cv 5? 


The H-characteristic and H-CUSP schemes introduce the modified Jacobian matrix 

dF 


Ah = 


dw h 


where Wh is the modified state vector 


^ P 
pu 

w h = pv 

pw 

\p» > 

Using transformations of the same kind a s those that were used for the standard Jacobian matrix, this can be 
decomposed as 

Ah ~ MhAhM^ 1 
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The eigenvalues of Ah, and also those of Ah , since they can be derived from Ah by a similarity transformation, 
are Q ,Q ,Q ,A + , and A” where 


A* = 



Thus A + and A change sign when Q = ±c5, and A + has the same sign as Q + cS, while A has the same 
sign as Q — c5. It is convenient to set 
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Also 


R h = DN h , R-^N^D- 1 

where D is the same diagonal scaling matrix as before, and 


/ 


N h = 


n x 

0 

n 2 

~ n y 

0 


n v 

-n 2 

0 

Hr 

0 


1 


n v 

~n x 

0 

0 


or n r a n x 
a + n 0 OL~n y 
a“n 2 

_ 2 

a 


a + n z 


a n 


with the inverse 


N, 


-1 


/ 

n y 

n z 

0 

V 0 


-27 an* 

—n z - 2y an x n y 
n y - 2yan x n z 

-cr~ . . 

o* n 

2/ia- 


n z - 2y<rn y n x 
-270-fly 

-n x - 2yan v n z 
:n v 


2ila^ n y 

-V + 

— n y 


2^ia~ 


— n y - 27 <rn z n a: 
n x — 27 <rn z n y 
- 27 <rn^ 




2/ia+ 
-1 


2fia“ / 


The formulas for the standard Jacobian matrix are recovered by setting a + = 1, a = — l,/i = ^(a + — a”) = 1, 
a = — |(a + -f a“) = 0. Correspondingly, Ah can now be represented as 

Ah = RhAhRh 1 

where the right eigenvectors of Ah are the columns of 

Rh — MhRh = MhDNh 


and the left eigenvectors of Ah are the rows of 

Ri l =ki l MZ l = N; l D- l MZ l 

These decompositions of A and Ah express every element in terms of velocities and metric quantities, while 
the density is completely eliminated from the formulas. 
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2a: RAE-2822 Airfoil 2b: NACA-0012 Airfoil 

Figure 2: O-Topology Meshes, 160x32 
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3a: C p after 25 Cycles. Wo,k 

Ci — 1.1227, Cd = 0.0460. 3b: Convergence. 


Figure 3: RAE-2822 Airfoil at Mach 0.750 and a = 3.0° 
H-characteristic Scheme. 



4a: C p after 25 Cycles. Wofk 

Ci = 1.1312, Cd = 0.0469. 4b: Convergence. 


Figure 4: RAE-2822 Airfoil at Mach 0.750 and a = 3.0° 
H-CUSP Scheme. 


5a: C p after 75 Cycles. 

Ci = 0.3620, Cd = 0.0230. 5b: Convergence. 


Figure 5: NACA-0012 Airfoil at Mach 0.800 and a = 1.25° 
H-characteristic Scheme. 



6a: C p after 35 Cycles. Wott 

Ci = 0.3654, Cd = 0.0232. 6b: Convergence 


Figure 6: NACA-0012 Airfoil at Mach 0.800 and a = 1.25° 
H-CUSP Scheme. 





7a: C p after 75 Cycles. 
Ct = 0.3818, C d = 0.0580. 



7b: Convergence. 


Figure 7: NACA-0012 Airfoil at Mach 0.850 and a = 1.0° 
H-characteristic Scheme. 



8a: C p after 35 Cycles. 
Ci = 0.3861, = 0.0582. 
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8b: Convergence. 


Figure 8: NACA-0012 Airfoil at Mach 0.850 and a = 1.0° 
H-CUSP Scheme. 




9c: 50.00% Span. 9d: 68.75% Span. 

Ci = 0.3262, C d = 0.0089. C, = 0.3195, C d = 0.0026. 


Figure 9: Onera M6 Wing. 

Mach 0.840, Angle of Attack 3.06°, 192x32x48 Mesh. 
C L = 0.3041, C D = 0.0131. 

H-CUSP scheme. 
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